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Tomonaga-Luttinger (T-L) theory predicts collective plasmon resonances in 1-D nanostructure 
conductors of finite length, that vary roughly in inverse proportion to the length of the structure. 
In-depth quantitative understanding of such resonances which have not been clearly identified in 
experiments so far, would be invaluable for future generations of nano-photonic and nano-electronic 
devices that employ 1-D conductors. Here we provide evidence of the plasmon resonances in a num¬ 
ber of representative 1-D finite carbon-based nanostructures using first-principle computational elec¬ 
tronic spectroscopy studies. Our special purpose real-space/real-time all-electron Time-Dependent 
Density-Functional Theory (TDDFT) simulator can perform excited-states calculations to obtain 
correct frequencies for known optical transitions, and capture various nanoscopic effects including 
collective plasmon excitations. The presence of 1-D plasmons is universally predicted by the various 
numerical experiments, which also demonstrate a phenomenon of resonance splitting. For the metal¬ 
lic carbon nanotubes under study, the plasmons are expected to be related to the T-L plasmons of 
infinitely long 1-D structures. 
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I. INTRODUCTION 


Nanoplasmonics is a field that has grown rapidly in 
the last few years 0 , 0 - Covering a range from terahertz 
through infrared to visible light frequencies, this field 
already offers numerous applications to electronics and 
photonics. In the visible and Near IR (NIR) range, nano¬ 
antennas (sj and nanoparticles have provided drastically 
enhanced coupling to electromagnetic waves [!]. All of 
this research work has taken advantage of plasmonic sur¬ 
face modes - either on nanoparticles, or on metallic sheets 
- and has led to some initial practical applications. How¬ 
ever, there are still significant fundamental limitations 
that must be overcome for further progress to be possi¬ 
ble in this field Q. While there is an extensive literature 
on 1-D plasmons, these have not been emphasized for ap¬ 
plications due to the difficulty of preparing devices and 
performing experiments. Yet, 1-D conductors are of in¬ 
terest since they can potentially exhibit lower losses than 
the 2-D ones. The chief signature of 1-D plasmons is a 
high-frequency excitation (i.e. energy resonance) that 
depends inversely on the length of the conductor. Exam¬ 
ples of nano-structures that can support 1-D plasmons 
are graphene nanoribbons (GNRs), single wall carbon 
nanotubes (SWCNTs) and linear carbon chains. Due to 
the very strong coulomb interactions in such 1-D systems, 
the low-energy excitations of interacting electrons cannot 
be adequately described by Landau’s Fermi liquid theory 
0, which has to be replaced by Tomonaga-Luttinger (T- 
L) liquid theory [7Hlu|. Specifically, it was shown that 
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metallic carbon nanotubes can be well described with 
the T-L theory [ll], [jj|■ A bosonized representation of 
the Hamiltonian led to the following conclusions: (i) the 
density of single particle states is a power function of the 
energy and vanishes at the Fermi energy; (ii) consistent 
with (i) the conductance of electrons tunneling into the 
CNT has a universal power law dependence on voltage 
and temperature; (iii) the electron states are separated 
into charge states and spin states [13]. In the simula¬ 
tions that follow spin/charge separation effects are ne¬ 
glected; (iv) collective boson states form charge density 
waves in the CNT (“T-L plasmons”) which propagate 
with a velocity vp = vp/g where vf is the Fermi ve¬ 
locity re 1 x 10 6 m/s, M and g is a parameter that de¬ 
pends on the strength of the Coulomb interaction and 
on screening by the electrostatic environment (the latter 
being a weak effect) [l2j . The predictions (i) and (ii) 
have been well verified by transport EMU and optical 
0 measurements. The parameter g is < 1 for T-L liq¬ 
uids and = 1 for Fermi liquids, vp is thus higher than 
Vp. An average plasmon velocity of the predicted mag¬ 
nitude was recently deduced from absorption measure¬ 
ments on enriched SWCNTs films |l8j. Ultra-fast time- 
dependent measurements evidenced ballistic electron os¬ 
cillations in isolated SWCNTs, but the propagation ve¬ 
locity was equal to the Fermi velocity, consistent with 
single particle excitations, not plasmons M- Detailed 
experimental confirmation of the predicted plasmon res¬ 
onance frequencies in isolated SWCNTs is thus still an 
unsolved problem. Estimates of g for CNTs yield values 
in the range 0.26 — 0.33 fa resulting in vp being in the 
range 3 — 4 x 10 6 m/s. Electromagnetic simulations have 
yielded vp up to 6.2 x 10 6 m/s [20j. Recently, T-L liquid 
behavior with g = 0.53 was verified in atomic gold chains 
by tunneling and optical measurements [2li ] . 
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The focus of this paper is the dynamic behavior of the 
collective medium in 1-D conductors and its resonances. 
A physical picture of the resonances can be gained by 
realizing the equivalence of the finite length 1-D conduc¬ 
tor to a transmission line [22|, [Hj]. The T-L plasmon 
waves resonate as they are reflected at the ends of the 
conductor and the fundamental resonance frequency will 
be given by / = v p /(2L) [24|] where L is the length of 
the conductor. Single quasi-particles similarly resonate 
at f = Vp/(2L). In our simulations we explore both 
of these types of resonances as well others that will be 
discussed in the text. 

First-principle calculations are becoming critical to 
provide evidence of plasmon resonances in 1-D carbon- 
based nanostructures. Such quantitative simulations are 
known to be challenging, since they should be both ca¬ 
pable to provide reliable information on the many-body 
excited states (beyond ground state theory), and to ad¬ 
dress large-scale computational needs of finite dimen¬ 
sional systems (beyond the solid-state unit-cell). The 
time dependent density functional theory (TDDFT) (25} 
alongside with the simple adiabatic local density approxi¬ 
mation (ALDA) for the many-body exchange-correlation 
term, has been very successful for providing accurate ab¬ 
sorption spectra of a large number of complex molecu¬ 
lar systems f2(J]. The real-time TDDFT approach in¬ 
troduced by Yabana and Bertsch HU, in particu¬ 
lar, combines both potential for parallel computing seal- 
ability, and an intuitive treatment of the real-time spec¬ 
troscopy that can deal with any form of excitations. Our 
in-house simulator, named NESSIE, is performing both 
all-electron real-space DFT ground state and real-time 
TDDFT excited states calculations. NESSIE benefits 
from the h igh-per formance capabilities of the FEAST 
eigensolver [29l 1301 ], which has also allowed to efficiently 
redesign various stages of the electronic structure numer¬ 
ical modeling process [3ll - l33l ]. Detailed information on 
our numerical real-space and real-time modeling frame¬ 
work is provided within the supporting document of this 
article. 

Two types of real-time TDDFT simulations are 
considered in our experiments. First, the spectroscopic 
information is obtained from linear response calculations 
after a short and weak polarized impulse is applied 
along the longitudinal or perpendicular direction of 
the 1-D nanostructure. Then, once resonances of 
interest are identified in the absorption spectrum, new 
time-dependent calculations are performed in response 
to a realistic stimulus, such as a laser tuned to each 
resonance frequency (e.g. sinusoidal excitation along 
the longitudinal direction). Such simulations aim at 
providing more detail on the 3-D electron dynamics 
of the particular resonances with relevant information 
about their nature. From atom and benzene-like chain 
structures to short carbon nanotubes, our numerical 
experiments discussed below progressively account for an 
increase in the physical and computational complexities 
of the 1-D atomic structure. 


II. CARBYNE 

Carbyne is a chain of carbon atoms that comprises ei¬ 
ther double or alternating single and triple atomic bonds. 
With recent progress in synthesis [34j , first-principle the¬ 
oretical investigations have become increasingly more rel¬ 
evant to study the various attractive physical properties 
of this material ]35s|. Previous studies of carbon chains 
using the TDDFT approach have been reported in Refs. 
13 (MI 81 . Our numerical simulations of NC 2 n N shown in 
Figure [Q are in good agreement with the experimental 
data found in Ref. [39|, where the strong longitudinal 
resonance shifts progressively towards the left of the spec¬ 
trum with longer chains (i.e. red shift). Their oscillator 
strength increases with the length of the chain indicat¬ 
ing that more electrons can participate in the plasmonlike 
collective excitations. Furthermore, the associated veloc¬ 
ity of the resonances is in the typical range o£ collective 
plasmons (as discussed further). References [3f| 38] also 
identify these excitations as collective plasmons. We note 
that the carbynes are not metallic and thus are not ex¬ 
pected to conform with the T-L theory. Nevertheless, 
carbynes show strong 1-D Plasmon excitations. More¬ 
over, our results show that these longitudinal resonances 
do not appear with a perpendicular excitation (see the 
NCgN plot at the top of Figure [l]). In contrast, one can 
identify other small resonances for instance, at ll.OleV, 
and at 14.68eV which do not depend on the length of the 
chain. These particular energy transitions are compara¬ 
tively close to the ionization potentials of the local atoms 
C and N, respectively at 11.26eV and 14.53eV [io| . 

The longitudinal resonances can be observed in more 
detail in Figure [2] (left plot) for the NCieN, HCieH, 
and HC 24 H chains. We note that the spectra produced 
for the Ci 6 chain using the H or N termination are 
very similar. Our results are also quantitatively close 
to the ones obtained in Ref. (38] for the HC 2 n H chain 
that makes use of the projector augmented-wave (PAW) 
pseudopotential and performs TDDFT linear response 
simulations in the frequency domain. The lowest energy 
plasmon resonances found in these results are also in 
agreement with earlier real-time TDDFT simulations 
performed in Ref. [37J using a smooth pseudopoten¬ 
tial. In comparison with the latter however, both the 
PAW and our full-core potential treatments reveal 
additional plasmonlike resonances at higher energies. 
These oscillations can be observed in the experimental 
data of the NC n N chains [39]. In our simulations, the 
splitting between these resonances becomes narrower for 
the HC 24 H and only two main peaks can be observed 
instead of three for HCigH. The isosurface snapshot 
plots in Figure [2] represent three additional real-time 
TDDFT simulation results that can help provide more 
insights on the time-dependent electron dynamics of 
the plasmonic resonances in HC 16 H. The strongest 
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FIG. 1: Computed absorption spectra of NC 211 N with lengths 
n = 1,2,4, 8 (using the single-triple bond alternation struc¬ 
ture and the optimized geometries reported in Ref. Hi)- The 
plot on the top for n = 4 is associated with the perpendicular 
response of a weak impulse applied along the same direction. 
All the other plots consider the response of an excitation ap¬ 
plied along the longitudinal direction of the chain. The figure 
insets illustrate the variation of the induced dipole obtained 
from our real-time TDDFT calculations and that are used to 
derive the absorption spectra of the corresponding structures. 
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3rd: E=4.41eV; f=I066.33 THz 


FIG. 2: The left plot includes a more detailed representation 
of the main resonances of NC 16 N observed in Fig. [l] it is 
also plotted alongside with the absorption spectra of the H- 
terminated carbon chains HC 24 H and HC 16 H (using the op¬ 
timized geometries reported in Ref. SI)- The latter counts 
three longitudinal resonances at 3.4eV, 3.94eV and 4.41eV. 
After new excitations of the structure at each one of these 
specific energy transitions (i.e. sinusoidal excitation at spe¬ 
cific frequencies), the figures on the right represent the 4D 
isosurface snapshots of the charge oscillation (i.e. the devia¬ 
tion of the charge density from the ground state - positive and 
negative deviations respectively in red and blue) all taken at a 
different time of the simulation where the dipole goes through 
a maximum and a minimum, respectively. 


resonance at the lowest energy transition gives rise to 
a distinguishable plasmon that can oscillate back and 
forth along the longitudinal direction. This collective 
particle excitation likely takes place within the first Van 
Hove singularity of the 1-D structure where the electrons 
are strongly coupled. It has the signature of a plasmon 
with a calculated velocity of v p i = 3.18 x 10 6 m/s which 
increases to v p i = 3.64 x 10 6 ? 7 i/s for HC24H. The 
results for the other two snapshot plots confirm the 
presence of additional lon gitu dinal resonance modes that 
were mentioned in Ref. [381 ]. with cosine and sine like 
envelope features for the plasmonic excitations. These 
two excitations have a different character in that the 
charge density alternates from positive to negative for 
each atom pair, analogous to optical phonon modes. 
Finally, we note that no single particle excitation at the 
Fermi velocity is observed since the carbon chain family 
considered here is semiconducting. 


III. ACENES AND POLY(P-PHENYLENE) 

Acenes and Poly(p-phenylene) (PPP) are polycyclic 
aromatic hydrocarbons consisting of linear benzene rings 
that are respectively fused or attached. They can also 
be thought of as the narrowest graphene nanoribbons 
(GNR) of finite lengths, associated with the zig-zag con¬ 
figuration for acenes (2-ZGNR) and the armchair one for 
PPP (3-AGNR). Understanding the properties of long 
acenes in particular, is the subject of active research [43| . 

Results for anthracene, tetracene and pentacene re¬ 
ported in Figure [3] confirm that the strong longitudinal 
UV resonances along the z-direction that appear respec¬ 
tively at 4.85eV, 4.30eV, 3.88eV, are in remarkable agree¬ 
ment with the experimental data in Ref. [45|. While 
TDDFT using traditional exchange-correlation function¬ 
als such as ALDA fail to describe the important low-lying 
excited states in acene compounds [H|, it appears capa¬ 
ble to accurately capture the main nanoscopic effects. 
These resonances are plasmonlike, since we note a pro¬ 
gressive shift towards the left of the spectrum with longer 
chains accompanied by an increase in oscillator strength. 
They are also absent from the results of excitations along 
the x and y directions. For the latter, in turn, one can 
identify two other resonances that stay independent of 
the chain length, at 11.22eV using an x-polarized impulse 
that equally excites all the carbon atoms in the y-z plane 
(energy relatively close to the ionization potential of C), 
and at 7.55eV using an y-polarized impulse that equally 
excites all the individual fused benzene rings (energy rel¬ 
atively close to the 7 r —> tt* transition in benzene). For 
the case of pentacene, in addition to a strong plasmon 
resonance at 3.88eV, we note a couple of weaker longitu¬ 
dinal resonances at 4.84eV and 5.25eV (see top right plot 
of Figure 0. The latter comes close to the experimen¬ 
tal HOMO-LUMO gap at 5.2eV [4(| (this gap is known 
to be severely underestimated by DFT at ~ leV). The 
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FIG. 3: The left plots present the computed absorption spec¬ 
tra of anthracene (N=3), tetracene (N=4), pentacene (N=5) 
associated with the responses of weak impulses separately ap¬ 
plied along the three main directions (regular bond lengths are 
used for the geometries: lc-c = 1.42 A and Zc-h = 1.09A). 
The top right plot focuses on the plasmonic response of an ex¬ 
citation along the longitudinal direction for pentacene (N=5), 
heptacene (N=7) and nonacene (N=9). The heptacene, in 
particular, counts two main resonances at 3.38eV and 4.41eV. 
The 4D isosurface snapshots represent the charge oscillations 
associated with these two energy transitions obtained at dif¬ 
ferent simulation times where the dipole response reaches a 
maximum and a minimum. 


results on heptacene and nonacene in Figure [3j clearly 
show that as the length of the structure increases, a sec¬ 
ond strong resonance appears. The isosurface snapshots 
representing the electron dynamics for the two main res¬ 
onances of heptacene confirm the presence of two plas- 
mons. In contrast to the case of carbon chains, this sec¬ 
ond resonance does not result from an additional longi¬ 
tudinal confinement mode, and the results indicate the 
presence of another channel (i.e. two confinement modes 
are then present in the transverse x 7 y plane). We note 
that similar double plasmon resonances have also been re¬ 
cently observed using semi-empirical simulations applied 
to GNR structures with variable widths [47j]. The cal¬ 
culated plasmon velocities associated with the first and 
second resonances of heptacene, v p i x = 2.77 x 10 6 m/s and 
Dpi 2 = 3.61 x 10 6 m/s, increase to v p i 1 = 3.18 x 10 6 m/s 
and v p i 2 = 4.08 x 10 6 m/s for nonacene. 

The case of finite PPP is shown in Figure @] using 
lengths of 5 and 10 unit cells. Excitations along the 
perpendicular directions (including the x-direction not 
shown here) provides results comparable to the ones 
obtained for acenes in Figure [5] However, we note one 
additional resonance at around 6.55eV for an excitation 
along the y perpendicular direction (curve denoted 
’5T’) which is also present (and amplified) with an 
excitation along the longitudinal z direction for both 


FIG. 4: The left plot presents the computed absorption spec¬ 
tra of a finite PPP with 5 and 10 unit cells (using the geome¬ 
tries reported in Ref. 0), obtained after excitations along 
the z longitudinal or y perpendicular direction. The result for 
the latter is noted 5T (for the 5 unit cells). The 4D isosurface 
snapshots on the right represent the charge oscillations asso¬ 
ciated with three main selected longitudinal resonances of the 
5 unit cell PPP. 

5 and 10-PPP. All the other longitudinal resonances 
at lower energies are absent from the ’5T’ curve. The 
isosurface snapshots of the charge oscillations, provide 
some useful information on the nature of the three main 
selected longitudinal resonances for 5-PPP. The strong 
first peak can be associated with a plasmon, while the 
second and third peaks reveal different characteristics 
related to “band-to-band” transitions. In particular, 
the second peak shifts from 3.86eV to 3.3eV for the 
long 10-PPP which is in very good agreement with the 
strong experimental absorption peak (related to the 
bandgap) found at 3.4eV for PPP 0. The position of 
the third peak, in turn, is not sensitive to variation in 
length. A closer look at the isosurface plots show that 
the charge oscillations at 6.55eV are somehow localized 
within each attached benzene ring. We note that a 
broader “band-to-band” transition region appears at 
around 5eV for both 5 and 10-PPP. Finally, the main 
plasmon resonance at 2.67eV i.e. v p i = 2.44 x 10 6 m/s 
shifts towards the left of the spectrum of the 10-PPP 
and, similarly to the case of long acene, splits into two 
peaks. The calculated plasmon velocities associated 
with the two resonances are v p i x = 3.97 x 10 6 to/s and 
Dpi 2 = 4.83 x 10 6 to/s. 


IV. SINGLE-WALL CARBON NANOTUBE 
(SWCNT) 

Figure [5] shows our results for simulations of a finite 
(3,3) armchair SWCNT. To stabilize the structure one 
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FIG. 5: The computed absorption spectra of the finite (3,3) 
armchair SWCNT with three different length 3, 5 and 7 unit 
cells (the C-C bond length is fixed at 1.44A and one addi¬ 
tional carbon ring of type ’A’ is added to the structures e.g. 
’H-AB-AB-AB-A-H’ for the 3 unit cells with Hydrogen termi¬ 
nation). The ’5T’ curve indicates the result from an excita¬ 
tion in the perpendicular direction, all other results consider 
an excitation along the longitudinal direction of tube. The 
4D isosurface snapshots on the right and below the spectrum 
represent the charge oscillations associated with five main se¬ 
lected longitudinal resonances of the 5 unit cell tube. They 
have all been taken at a different time of the simulation where 
the dipole goes through a maximum, and the figures include 
both side and edge views of the nanotube. 


additional carbon ring, and hydrogen atoms were added 
at the ends. The H atoms are expected to have only mi¬ 
nor effects on the resonances observed. Armchair SWC- 
NTs are known to exhibit metallic conduction [ldj and, 
as mentioned in the introduction, they can be described 
by the T-L theory in the infinite length limit [Tl], [l2| ■ 
We simulated tubes of three different lengths, 3, 5, and 
7 “unit cells”, respectively. The simulations were per¬ 
formed for an E-field parallel to the z-axis of the tubes 
(the long axis), marked ’3’, ’5’ and ’7’, as well as perpen¬ 
dicular to that axis (for the 5 unit cell case, marked ’5T’). 
We first distinguish resonances that do not change with 
the length (L) of the tube. The peak at about 4.44eV 
can be identified as being due to the 7r-plasmon (E n in 
the Figure). Note the characteristic distribution of the 
oscillating charges, following the bonds, in the snapshot 
displayed below the spectra. A broad band with a peak 
at ~ 15eV due to the 7r + cr plasmon can also be observed 
in the spectra (observations are included in the support¬ 
ing information document). Since the above plasmons 
are basically of the bulk (3-D) type, they can also be 
excited by a field in the direction perpendicular to the 
tube. Both of these resonances are well known for SWC- 
NTs and have been detected experimentally by Electron 
Energy Loss (EEL) spectroscopy [501 as well as by opti¬ 
cal spectroscopy (the 7r-plasmon) [5ll [521]. We note that 


the 7r-plasmon resonance for perpendicular excitation is 
shifted to a higher energy, 5.8eV, in qualitative agree¬ 
ment with the experimental data. The 7r-plasmon res¬ 
onance for the 3 unit cell case is not so easy to iden¬ 
tify, see further discussion below. There is also a reso¬ 
nance for all values of L at about 3.25eV (Ebb in Figure 
El. This resonance energy agrees with experimental op¬ 
tical data and simulations 15 ll 1531 ] and can be identified 
as a “band-to-band” resonance for infinitely long tubes 
and for higher energy bands. The simulations show that 
this peak does not vary with L, and also that it can not 
be excited perpendicular to the axis, as expected for a 
band-to band transition. The charge distribution under 
sinusoidal excitation shown in Figure El clearly lacks the 
collective character evident in the plasmon resonances. 

One type of resonance that does vary with L occurs in 
the lowest energy range, from about 1.4 to 2.2eV (E e - 
in the Figure). Calculating a velocity for this resonance 
we find a value 0.98 x 10 6 m/s and 1.18 x 10 6 m/s respec¬ 
tively for the 5 and 7 unit cell cases, which is close to the 
Fermi velocity. We note that this resonance can only be 
excited by a field parallel to the tube axis, as expected if 
it is due to electrons resonating from one end of the tube 
to the other. This is the type of single particle excita¬ 
tion resonance that was measured in Ref. [l9|. Further 
evidence for this interpretation is obtained by inspecting 
the top 4-D isosurface snapshots to the right of the spec¬ 
tra (’Single Excitation - E e - = 1.63eV’), which show the 
periodic oscillation of the charge density. 

Finally, we find two strong peaks at energies of 3.69 — 
3.89eV (a split peak for the longest tube) and 3.89eV 
(marked E p i i in the Figure) that vary with L as expected 
for T-L plasmon type resonances. These resonances can 
only be excited with a parallel field as was confirmed 
by a simulation of the 5 unit cell tube with perpendic¬ 
ular excitation (marked ’5T’). One of the series of 4-D 
snapshots (Collective Oscillation - E p i i = 3.89eV) shows 
how the charge density waves travel back and forth on 
the tube under sinusoidal excitation at the resonance fre¬ 
quency. There is a related plasmon resonance at a higher 
energy ( E p i 2 at 4.76eV) for the 5-unit cell case. Thus 
we find split plasmon resonances for both cases - the 5 
unit and 7 unit cell ones - with a much smaller split for 
the longer tube. The calculated velocities for the two 
T-L plasmons for the 5 unit cell case are 2.35 x 10 6 m/s 
and 2.87 x 10 6 m/s, and they increase to 3.12 x 10 6 m/s 
and 3.29 x 10 6 m/s for the the 7 unit cell. The side 
view snapshots of the charge distribution in Figure E] can 
now be used to discuss further the different characters of 
the resonances E e ~, E p q and E p i 2 . The charge density 
for the single (quasi-particle) excitation ( E e -) is spread 
across the cross-section. In contrast, the plasmon excita¬ 
tions concentrate the charges around the periphery and 
at the ends of the tube as expected for a boson-type mode 
The split between the two plasmon frequen¬ 
cies can be interpreted as related to two different modes 
quantized around the periphery. 

Given our interpretation of the peaks from 3.69eV to 
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FIG. 6: Summary of plasmon velocities computed for all the 
1-D carbon nanostructures using different lengths (numeri¬ 
cal values are provided in supporting information document). 
Two main plasmon modes are represented (except for the car- 
byne and the 5-PPP cases). We note the universal monotonic 
increase of the plasmon velocities for longer structures, while 
the less regular behavior in the 15-20A range (in particular 
for the 5-PPP and the 7-unit cell (3,3) SWCNT) is likely due 
to coupling effects between “band-to-band” transitions and 
the plasmons. 

3.89eV as T-L plasmon resonances, we would expect a 
similar T-L plasmon resonance for the 3 unit cell case at 
a somewhat higher energy, where we see two peaks. As 
indicated earlier, this brings us to the expected energy 
of the 7r-plasnron, however, and a possible interpretation 
of the two peaks at 4.32eV and 4.81eV for the 3 unit 
cell tube is that the T-L plasmon and the 7 r-plasmon are 
coupled. 

Conclusions 

This paper has simulated molecular absorption spectra 
employing accurate all-electron real-time TDDFT sim¬ 
ulations. The range in size has been extended from 
small molecules such as C 2 H 2 or benzene to carbon 
nano-structures that are equivalent to 1-D conductors 
with finite lengths. Excellent agreement with experi¬ 
mental spectra is obtained when such data is available. 


We find universal features for all structures investigated 
which include carbon chains, narrow armchair and zigzag 
graphene nanoribbons (i.e. acenes and PPP), as well as 
short carbon nanotubes. All structures show collective 
plasmon oscillations characterized by resonant frequen¬ 
cies that vary roughly inversely with the length of the 
structure. For the metallic structures (in our case only 
the SWCNTs) the plasmons are expected to be related 
to the T-L plasmons of infinitely long 1-D structures. 
Another notable characteristic is that the main plasmon 
resonance is split into two components due to transverse 
quantum confinement effects (except for the case of car- 
byne that presents two to three longitudinal resonance 
modes). The simulations allow vivid 4-D visualizations of 
the electron charge distributions for different types of res¬ 
onant modes, which then provide more insights on their 
nature. As summarized in Figure | 6 l the velocity of plas¬ 
mon propagation shows a monotonic increase with the 
length of the structure and it is expected that the plas¬ 
mon velocity will asymptotically become independent of 
length. Further optimization of our present high perfor¬ 
mance computing techniques should allow the extension 
of the simulations to longer structures (up to tens of unit 
cells for SWCNT), and lead to accurate predicted data 
that will guide future photonic and electronic applica¬ 
tions of nanostructures over a wide frequency range, from 
visible to terahertz. 


Supplementary information 

The supplementary document to this article includes 
detailed information about the theoretical models and 
numerical methods used in this study, additional in¬ 
formation data about the 7 r and n + a resonances in 
SWCNT, and a table of numerical values for the plas¬ 
mon velocities. 
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1 Methods - Summary 


From physics to algorithms, an account of the multi-step modeling process used by our in-house 
real-space and real-time NESSIE simulator is summarized in the following. 

1.1 Ground-state calculations. 

Our physical model to calculate the ground-state wave functions is at the level of Density Func¬ 
tional Theory (DFT) associated with the Kohn-Sham equations. For this study, the many body 
exchange-correlation term is defined using a standard local density approximation (LDA) formula 
[1J. The electronic Hamiltonian accounts for the full core potential (also called all-electron cal¬ 
culations). In contrast to pseudopotential methods in particular, this approach provides higher 
accuracy, consistency, and higher locality in the discretization stage with better potentiality for 
parallelism. The number of electrons in our simulations range from a few tens for small molecules 
(e.g. 42 electrons for benzene) to a few hundreds for most of the nanostructures considered here 
(e.g. 552 electrons for the (3,3) SWCNT with 7 unit cells). The finite element method is used to 
discretize the Hamiltonian system using tetrahedral elements with either quadratic P2 or cubic P3 
basis functions (P2 is preferred for our TDDFT linear response simulations using weak perturbation 
- it is computationally less expensive and essentially provides the same results). The real-space 
simulation domain is taken large enough to minimize any unwanted scattering from artificial bound¬ 
aries. More details on our FEM real-space mesh approach are presented in the next Section. In 
turn, the significant challenges posed by the NESSIE’s large-sparse matrix computations result¬ 
ing from the real-space mesh discretization have been addressed recently with the development of 
our FEAST eigensolver 121 El- FEAST has considerably broadened the perspectives for enabling 
high-performance parallel large-scale calculations. In particular, we operate FEAST beyond the 
current “black-box” solver by considering a fundamental redesign of the electronic structure nu¬ 
merical modeling, including: (i) a reformulation of the real-space muffin-tin domain-decomposition 
problem for solving the first-principle all-electron problem exactly [4], (ii) a fundamental non-linear 
FEAST numerical technique for addressing the self-consistent (SCF) DFT ground-state problem 
that represents a robust alternative to traditional SCF methods using iterative procedure with 
density mixing schemes [5]. 

1.2 Excited-state calculations. 

Within the real-time TDDFT framework, all the occupied N e ground-state wave functions {<pi, ..., 4>n b } 
are propagated (non-linearly) in time by solving a time-dependent Kohn-Sham type equation. The 
absorption spectrum is obtained after a short polarized impulse is applied in any given direction r/ 
(e.g. x,y or z ) of the molecular system i.e. = 0 + ) = ex.p(—ikr])(j)j(r) [6j. Using the solution 

of the electron density at each time step i.e. n(r, t) = 2 J^/=i |Vh( r >^)| 2 ) the induced dipole can 
then be computed in the time-domain d{t) = f drn(r, t)r /, and the imaginary part of its Fourier 
transform provides the dipole strength function and the absorption spectrum. Examples on the 
dipole variations and more details on the calculations of the dipole strength function are presented 
in the next Section. The real-time TDDFT technique is often used to address linear responses 
(e.g. using weak momentum k ), but the framework can also generally address any forms of non¬ 
linear responses. In our simulations in particular, a longitudinal sinusoidal excitation is often used 
to observe the dynamics of the charge oscillations at a given resonance frequency. The real-time 
propagation is performed using either a conventional Crank-Nicolson (CN) scheme or a spectral 
decomposition technique where FEAST is used to diagonalize the time-dependent Hamiltonian at 
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each time step fTj. This FEAST spectral decomposition technique is well-suited for parallel com¬ 
puting, and it allows to consider larger time-steps in our simulations A t = 20 as (vs A t = 2 as using 
CN). The total simulation time that has been used varies between lOfs for the carbynes to 15fs for 
the SWCNTs. 

2 Methods - Practical considerations with applications to benzene 

Supplementary information on the modeling procedure used in NESSIE is provided here. We use 
the example of the benzene molecule. 

2.1 Real-space discretization 

As specified by an input file, our simulator NESSIE reads the atom coordinates, creates a 3D 
mesh using tetrahedra, and it builds the FEM Hamiltonian and mass (i.e. basis function overlap) 
matrices using quadratic P2 or cubic P3 precision. 

The 3D FEM mesh is constructed using a muffin-tin decomposition technique that involves two 
steps: (i) a 3D atom-centered mesh which is highly refined around the nucleus to capture the full 
core potential, and (ii) a much coarser 3D interstitial mesh that connects all the atom-centered holes 
and extends to the computational domain boundary. Figure ISll illustrates the essence of our FEM 
muffin-tin domain decomposition for the benzene molecule. The size of the full 3D (reconstructed) 



Figure SI: Using a muffin-tin domain-decomposition method, the whole simulation domain P is 
separated into multiple atom-centered regions (i.e. muffins on the right) and one large interstitial 
region (on the left). The figure in the middle represents a 2D section of local finite element dis¬ 
cretization using a coarse interstitial mesh (represented only partially here) connecting all of the 
atoms of a benzene molecule. 

sparse system matrix for benzene is ~ 60 K using P2, and ~ 140 K using P3. In TDDFT calculations 
using a weak perturbation, a P2 FEM basis provides results similar to P3. Table ISll reports the 
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Carbyne 

Acenes 

PPP 

SWCNT 


HCi 6 H 

HC 24 H 

Penta. 

Hepta. 

Nona. 

5PPP 

10PPP 

5-unit 

7-unit 

# electrons 

98 

146 

146 

198 

250 

202 

402 

408 

552 

Size H matrix 

65,357 

93,115 

138,926 

180,272 

223,813 

191,643 

368,133 

302,282 

393,410 


Table SI: Hamiltonian system size P2 FEM discretization by NESSIE for various 1-D carbon 
nanostructures. The system matrices are very sparse with an average of 15 to 20 non-zero elements 
per row. 


size of the resulting P2 FEM Hamiltonian system matrices and the total number of electrons for 
most of the 1-D carbon nanostructures we have considered in our simulations. 

The real-space simulation domain H is taken large enough to minimize any unwanted scattering 
from artificial boundaries. The wave functions are set equal to zero at the boundaries (i.e. Dirich- 
let boundary conditions for Kohn-Sham), and radiative boundary conditions on the electrostatics 
potential are used for solving the Poisson equation. The coupled DFT-Kohn-Sham/LDA/Poisson 
problem is solved self-consistently using the non-linear version of the FEAST solver. In particu¬ 
lar, the lowest N e /2 states are computed (N e denoting the number of electrons, and the factor 2 
accounts for the spin). 

2.2 Real-time calculations and absorption spectra 

Two perturbations of the ground state solutions are considered in the simulation (i) a weak short 
impulse along a given direction (longitudinal or perpendicular to the structure), (ii) a weak si¬ 
nusoidal stimulus at a given resonance frequency. In both cases, the ground-state occupied wave 
functions (i.e. electrons) need to evolve in time. At each new time step t + At, the Poisson equation 
is solved, and the time-dependent ALDA Hamiltonian is updated. We are also making use of a 
simple corrector-predictor schemes to evaluate the Hamiltonian at t +Aj/2. For the real-time prop¬ 
agation, two strategies have been used: (i) a Crank-Nicolson (CN) scheme that requires solving a 
linear system at each time step (number of right-hand-side is equal to N e / 2); (ii) a spectral decom¬ 
position scheme that requires the partial diagonalization of the large Hamiltonian at each time step. 
The latter direct diagonalization technique can become a viable alternative to CN (i.e. potentially 
capable of both higher-scalability and better accuracy) only if used within a parallel computing 
environment and by taking advantage of various intrinsic properties of the FEAST solver (e.g. the 
eigenvector subspace computed in the current time-step can be used as a good initial guess for 
computing the next one). 

Thereafter, at each time step we compute the 3d electron density (sum of amplitude square 
of the occupied wave functions that are propagated). The variation of the induced dipole is then 
computed relative to the dipole of the ground state and the center of mass tq = (xo, 2 /Oi^o) of the 
nanostructure. For example, the response along the longitudinal z direction is given by: 

Sd z (t) = f dr(n(r,t) - n(r,0))(z - z 0 ). 

Jn 

For the case of a short polarized impulse initial perturbation (using a weak momentum k), the 
dipole strength function S is obtained by taking the imaginary part of the Fourier transformation 
of the dipole moment, i.e. 


S(E) 


2 mE c 
nkfi 3 


dt 5d(t)e lEt / n w(t) 
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where w(t) represents a damping factor that introduces artificial dissipation for suppressing the 
noise in the Fourier transform (we choose w(t) = e _7i / r ). Typically in our simulations, we set 
7 = 3, and k = 0.01/L where L is the dimension of the computational domain along the polarization 
direction. For the example of the benzene molecule, Figure IS2l presents the results obtained for the 
dipole variations and their associated absorption spectrum, after an excitation is applied along the 
three main directions of the molecule. The mean of the three solutions is plotted in Figure [S3] and 
it can be quantitatively compared with the experimental data. 







Figure S2: Dipole variations and corresponding absorption spectrum after excitation along the three 
main direction of the benzene molecule. The CN scheme is used for real-time propagation with 
time-step A t = 5as and a total time of lOfs (2000 time steps). 
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Figure S3: Calculated absorption spectrum of the benzene molecule using NESSIE and comparison 
with experimental data. 
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3 Supplementary Figure for SWCNT 


Figure 5 shows the simulation results for a finite (3,3) armchair SWCNTs including the it + a 
bulk-type plasmon peak - broad band at ~ 15eV. We note that both it and it + a resonances for 
a perpendicular excitation are shifted to higher energy values in qualitative agreement with the 
experimental data (e.g. Kramberger C. et ah, Phys. Rev. Lett., 100, 119803, 2008). 



E(eV) 


Figure S4: Calculated absorption spectra of the finite (3,3) armchair SWCNT with three different 
length 3, 5 and 7 unit cells (the C-C bond length is fixed at l.ffA and one additional carbon ring 
of type ’A’ is added to the structures e.g. H-AB-AB-AB-A-H 1 for the 3 unit cells with Hydrogen 
termination). The ’5T’ curve indicates the result from an excitation in the perpendicular direction, 
all other results consider an excitation along the longitudinal direction of tube. 


4 Supplementary Table for Plasmon velocities 



Carbyne 

Acenes 

PPP 

SWCNT 

HCieH 

HC 24 H 

Penta. 

Hepta. 

Nona. 

5PPP 

10PPP 

5-unit 

7-unit 


19.34 

29.67 

12.29 

16.94 

21.78 

18.86 

39.81 

12.48 

17.48 

E Ph (eV) 

3.40 

2.53 

3.88 

3.38 

3.02 

2.67 

2.06 

3.89 

3.69 

v p i 1 (10 6 m/s) 

3.18 

3.64 

2.31 

2.77 

3.18 

2.44 

3.97 

2.35 

3.12 

E p i 2 (eV) 



4.84 

4.41 

3.87 


2.51 

4.76 

3.89 

v p i 9 (10 6 m/s) 



2.88 

3.61 

4.08 


4.83 

2.87 

3.29 


Table S2: Numerical values of T-L plasmon velocities computed for all the 1-D carbon nanostruc¬ 
tures using different lengths. Two main plasmon modes are reported (except for the carbyne and 
the 5PPP cases). 
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